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ABSTRACT 


A method is develoned for digital simulation of linear time- 
invariant dynamic systems with lumped parameters and time delays. 
Ordinarily, such systems can be described by a linear matrix differential- 
difference equation, which can be transformed to an infinite-dimensional 
difference equation whose solution is obtained in a recursive way. 


As the present method depends on the accuracy of evaluation of the 
matrix exponential, a simple computational, procedure based on the 
truncation of the infinite series for e is described. 

In addition, an algorithm is given that ensures that the 
transient state of an unforced linear time-invariant dynamic system with 
zero time delay is calculated to a specified accuracy. 


Several sample problems are included. 
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CHAPTER 1 


INTRODUCTION 


1-1 Description of the problem 


This report presents a method for the simulation of linear time- 
invariant dynamic systems with lumped parameters and time delays. 

In many industrial processes one often encounters a type of time 
delay called “transportation lag". This kind of delay is generated when 
process materials move from one point in a process to another point 
without any appreciable change taking place in the properties or 
characteristics of the process materials. Such delays may be caused by the 
flow of fluids through pipes, or by the motion of webs or filaments. 
Systems such as distillation columns and long heat exchangers are 
characterized by a multitude of small lags, which have an effect somewhat 
similar to that of time delays. The effects are not identical; however, 
some insight may be gained by using time delays models, The control of 
composition in a chemical reactor has been selected as a typical problem 
and this is depicted in section 5-2. 

Models having delays often arise in the study of systems with a 
mixture of lumped and distributed elements. An interesting form of 
topological representation suitable for such systems has been invented by 
Prof. H. M. Paynter at M.I.T., and is called the bond graph. Rosenberg (17) 
and Auslander q)* describe its use in modeling in some detail. 


Many other physical systems, such as electrical, mechanical and 


* Numbers in parenthesis refer to items in the bibliography. 


hydraulic transmission lines, and certain types of structural problems, 
are good examples of distributed systems which can be modeled using the 
delay operator. These systems often are analyzed as two-port chains, and 
usually the equations are slightly more involved than the type treated in 
this report. It is suggested that the reader interested in these kind of 
problems consult Koepcke (9) and Vaughn (20), as well as any standard text 


treating transmission phenomena, 
1-2 Formulation of the approach 


As an extension to the use of ordinary differential equations 
which arise when the future behavior of the system depends only upon its 
present state and not upon its past history, many systems that include 
time delays can be described by a linear matrix differential-difference 
equation. That is, the system is described by 


° n m 
X(t) = A, X(t - T,) + D, Ut -T,) , 
ee 1 bs j 


where X and U are the state and input vectors, respectively and T, and Ty 


are some fixed delay times, A, are a set of n x n matrices, and D, are a 


1 j 

set of n x r matrices, Techniques such as the direct method of lyapunov or 
laplace transforms can be used in the analysis of the equation. However, 
the use of these techniques frequently requires extensive computation, and 
for that reason they are not practical for hand analysis. At this step, 
designers and analysts are forced to rely on the digital computer as a 
computing aid, 


Because matrix manipulations are so convenient to implement on a 


digital computer, many existing dynamic systems programs are based on a 


matrix formulation of the problem, This convenience, together with the 
inherent elegance of the matrix approach, is helping to promote its 
acceptance among systems theorists. 

This report analyzes systems governed by the following differential- 
difference equation, for which it is desired to have a time sampled version 


of the state response: 


X(t) = A X(t) +B X(t = 7) + DUC) + D,We - TD 
where 

T = time delay, 

X(t) = (n x 1) vector. It is called the state vector, 

U(t) = (r x 1) vector. It is the forcing signal or 
input vector, and it is assumed to be constant 
between samples. 

A, B = (n x n) constant coefficient matrices, 


Die Dy = (n x r) constant driving matrices, 


Koepcke (9) shows that the equivalent difference equation is (see 
section 3-1) 


X(t +1) =) [0 (t) X(t = aNt) + A, (1) UCt - aNt)] 
i=0 


where N = 2 » and e, and dy are called plant transition matrices and 
control transition matrices respectively. 

The accuracy of evaluation of these sets of transition matrices 
depends upon the accuracy of evaluation of the matrix exponential, In 
section 2~3 a simple procedure based on the truncation of the infinite 
series of Pg (11,6), which guarantees a specified accuracy in the matrix 


exponential, is described. 


Also, a procedure is developed (21) to ensure that the calculated 
transient state of unforced linear time~invariant dynamic systems with 
zero time delay, is accurate to a specified tolerance. 

Several sample problems are presented to demonstrate the 


computation techniques. 


1-3 Application of results in dynamic simulation 


The two sets of simulators deduced throughout the development of 
this work, were tested on the time-shared IBM 7094 operated by Project MAC, 
and the entire operation, input and output, was carried out at an IBM 1050 
remote console typewriter. The algorithms will be part of the ENPORT 
Project which is being carried out at the mechanical engineering department 
under the direction of Professor Rosenberg, 

ENPORT is a digital computer program that accepts a bond graph 
description of a dynamic system and produces its time response. Work is 
being done on the theory of bond graphs, and a systematic graphical 
method has been developed for generating the state differential equations. 
ENPORT is organized in such a way that a broad class of nonlinear, active 
and passive, mixed energy~type systems can be handled. 

The wakelike nature of certain types of distributed systems make 
simulation by means of the digital computer, with its ability to exactly 
model the time delay operator, very natural. A simulation method based on 


delay-bond modeling has been developed by Auslander (1). 


CHAPTER 2 


DYNAMIC RESPONSE OF LINEAR TIME-INVARIANT SYSTEMS 


The analysis of many systems problems encountered in scientific 
and engineering investigations can be performed by either one of two 
major approaches. The essentially block diagram approach, involves the 
determination of the transfer characteristics of the system components 
and the overall transfer characteristics. The second approach is based 
primarily upon the characterization of a system by a number of coupled 
first order differential equations which govern the behavior of the state 
variables. This technique is often implemented with the aid of a state 


variable diagram and is referred to as the state-variable approach. 


2-1 System Characterization by State Variables 


From the point of view of system analysis it is convenient to 
classify the variables which characterize or are associated with any 
system into (1) input, or forcing signals, Ui, which in essence represent 
the stimuli generated by systems other than the one under investigation 
and which influence the system behavior; (2) output, or response, 
variables Yi, which describe those aspects of system behavior that are 
of interest to the investigator; and (3) state variables Xi, which 
characterize the dynamic behavior of the system under investigation, 

One way of defining state variables is by making use of the state 
variable diagram. A state variable diagram is made up of integrators, 
coefficients and summing devices. It describes the relationships among 


the state variables and provide physical interpretations of them. The 


outputs of the integrators denote the state variables. 

For continuous-time systems the state variable diagram is the 
same as the analog-computer simulation diagram. The state variable 
diagram may be derived from the overall transfer function of the system 
in three different ways (1) direct programming, (2) parallel programming, 
and (3) iterative programming. These methods are later ilustrated in the 
chapter corresponding to the solution to sample problems. Further 
information can be obtained from Tou (19), Schwarz and Friedland (18) and 


Ogata (15). 
2-2 Digital Solution of the Matrix Differential Equation 


A linear time-invariant system or process can be described by a set 
of first order linear differential equations with constant coefficients, 


which may be expressed in matrix form as 


X(t) = A X(t) + D U(t) (2.1) 
where 

A is the coefficient matrix 

Dis the driving matrix 


X is the state variable vector 


Ic: 


is the state forcing signal vector 


By analogy to the scalar case, the solution of.eq. (2.1) is 


T 
eA(Tt - T) 
t 
o 


A(T - t,) 


X(T) =e X(t.) +| D U(t) dt (2.2) 


with the initial conditions given by X(t,). 
For simplicity let ty = 0, and let us define 


AT 


o(T) =e (2.3) 


as the transition matrix of the process. An equivalent name is the matrix 


exponential. 


Therefore eq. (2.2) can be reduced to 
That 
X(T) = @(T) X(0) + (7) fe D U(t) dt (2.4) 
0 
If T is small compared to the shortest period of interest in U(t), 


U(t) may be approximated over the region by U(0). 


Then eq. (2.4) becomes 


X(T) = (T) X(0) + on jew ar] D u(d) (2.5) 
By integration of the series of e** 
fem dt = ati ~ o(-T)] (2.6) 
Thus 
1 


X(T) = 0(T) X(0) + O(T) A [1 - #(-T)] D UO) (2.7) 


Let us define 


nl _ AT ,-1 (nat 


A(t) = [eA? 4 }D (2.9) 


as the control transition matrix. 


AT 


From the series definition of e, it is observed that 


at eAT = e Al at 


Therefore, eq. (2.9) becomes 


Act « teAT grt. @AT Q-AT gol) 
A(t) = feA? am? ~ am) 


or 


act) = ((e“? - 1) av4] (2.10) 


Thus eq. (2.8) can finally be written as 


xcr) = eA™ xo) + (eA? - 2) a4] dD UCo) , (2.11) 
or 
X(T) = 0(7) X(0) + A(T) U(0) (2.12) 


In general eq. (2.12) can be expressed as 


X(KF1T) = O(T) X(KT) + A(T) UT) (2.13) 
which indicates that the state vector of the process after a particular 
interval depends upon the previous vector and also depends upon the 
forcing vector evaluated at the previous time. 

There are several methods available for computing the closed form 
expression for or. either as a special case of the study of the functions 
of a matrix or by a purely algebraic mathod based on the Laplace Transform. 
It is suggested, for those interested in these schemes, that they consult 


Ogata (15), Zadeh and Desoer (23), or Bellman (2). 


2-3 Digital Evaluation of the Matrix Exponential 


eft is given by 
2 
AT B B, B B, B 
e me el+B+o( 5) +3057) +o (2.14) 


note that each term in parenthesis is equal to the previous term. This 
provides a convenient recursion scheme, 

To ensure a reasonable truncation of the series, it is necessary 
to judge the convergence of the series. The norm of a matrix A is a real, 
non-negative number denoted by llall » that gives a measure of the size of 
the matrix elements. 


Let 


AT OM +R 


e(t) #e 
where M is the truncated matrix which is an approximation of oe (see 


reference 11) 


K i 
w=) Goo (2,15) 
i=0 
and R is the remainder matrix 
oo A i 
R= J eu (2.16) 
i=K+1 


If each element in the matrix ad is required with an accuracy 


of "d" significant digits, then 


In ,ls 1079 |, 4| (2.17) 
where Thy and ™4 are elements of matrices R and M respectively. 
Let us define the norm of matrix A as: 
tall = min{max (la,,l1 max (hla, |} (2.18) 
ij 43 i 
For this norm, we have 
Va sllcllallllel (2.19) 
Ja, ,l<hal (2.20) 
and 
Vale lle Jal + [sll (2.21) 


Then, it follows that 
4 ax)" + Halt xt 
ley lik i! ley i! (2.22) 
i=K+1 i=K+1 


if the same norm is applied to the remainder matrix R. 


Upon expansion of eq. (2.22) 


10 


| yX+2 K+1 | jk+2 K+2 
| j < Alt + Alor + (2.23) 
= (K+1)! (K+2)! 


and, calling e the ratio of the second term to the first 


al X*2 Kt2 


_ cee2yt | alle (2.24) 
Jalk*2 hth K+2 ° 


(Dt 


Therefore 


als a (2.25) 


Making the substitution of eq. (2.19) into eq.(2.23), it follows 


that 
Ie, slog" Hal, Ilcary* {I (Ar)? + 
K! — (K+2) (K+1) 
, lank ll 3 
Rt Rao » (2.26) 
or 


ls NaI (Lal. Lac_ILac_I, 


K+2 K+1 


, Lac_Shar_I bac, ean ] « (2.27) 


K+3, -K4+2 0 K+1 


Thus 


Iz,,ls| Cale I Las, [ar [14 [ar_] 


K+2 
4! At il At | 
Kes KW --)| (2.28) 


Now, because any factor of the form lac for a>2 is always less 


than e€, by eq. (2.24), then 


li 


K 2 3 4 
Iry,ls Hoy" H lac {leetetete + -—-} (2.29) 


If e<l, eq. (2.29) takes the form 


Ir, [<lcax*t (la Eas: ee (2.30) 


K! K+l ° l-« 

This equation is suggested by Everling (6) as the upper bound in 
the remainder matrix R. 

In order to initialize the procedure, a certain K has to be 
chosen, but this K cannot be arbitrary, because it may happen that e>l, 


and relation (2.30) would not hold any more. 


This situation can be solved using eq. (2.25); thus 


xo WA 
="'e€ 


In order to ensure that e<1/2, the initial condition for K should 

be 
K> 2 Harl (2.31) 

However, it is possible that || Atl be less than 1/2; then K would 
be less than one. So, in order to avoid this possibility, an initial value 
of K can be obtained from 

K = max [ 2larl, 2 ] (2.32) 

At this point, Everling (6) suggests that K be incremented by 
half of its initial value, in the course of iteration. 

Although the matrix series approach for the evaluation of the 
transition matrix is suitable for digital computation, the disadventage 
stems from the convergence requirements for the series Pa so it would 
be desirable to speed the computation. 


This can be done recalling the basic relationship 


12 


a 
eft a [| (2.33) 
where a is chosen using the following expression 
a= 2s max(|a |" (2.34) 
= 1 
i,j 


where 8 is the smallest integer allowed. 


At/a 


The idea is to compute e » because the norm of At/a is smaller 
than At, and the series will converge faster. Once the addition of the 
corresponding elements in the matrix terms of the infinite series is done, 
all that is required is to raise the result to the power a. The last step 
involves very few matrix multiplications, because a is a power of 2; for 
example, if a = 32 only 5 matrix multiplications are performed at the end 
of the computation. 


The steps presented in this section are summarized in a flow 


diagram in chapter 4. 
2-4 Error bounds in the transient response 


Although the matrix et can be obtained within prescribed 
accuracy, the truncation error of the matrix series, and the roundoff 
error do propagate in the state vector with increasing time. 

It is desirable, therefore to derive recursion relations which 
bound the propagated error due to these sources. Whitney (21) suggests 
one method. 


The homogeneous case of eq. (2.13) is 
X(K+1 T) = 0(T) X(K T) (2.35) 


If eq. (2.15) is used in place of $(T), the numerical calculation 


13 


reads 

X, (K+1 T) =M X, (Kk T) (2.36) 
where X,(K+1 T) is the perturbed state vector obtained from numerical 
calculation. 


The propagated error at time (K+1l) T due to the approximate M 


is 
E(K+1 T) = X(K+1 T) - X, (K+ T) (2.37) 
Rewriting eq. (2.35) and substracting eq. (2.36) from it yields 
X(K+1 T) - X,(K+1 T) = [MARI[X,(K T) + E(K T)] 
-M X,(K T) » (2.38) 
or 
E(K+1 T) = [M+ RJE(K T) + R X,(K T) (2.39) 
From eq. (2.17) 
In,,l< 10° |m, | (2.17) 
We can define 
R, = In, ,! I (2.40) 


where I is a matrix each of whose elements is unity. Replacing R with R, 


in (2.39), we obtain the running error bound for E(K+l T), that is 


E(K+l T) = [M+ R] E(K T) + R, X,(K T) (2.41) 


The computation may be initialized assuming E(0) is zero. 
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CHAPTER 3 


DYNAMIC RESPONSE OF LINEAR TIME-INVARIANT SYSTEMS 


WITH LUMPED PARAMETERS AND TIME DELAYS 


It has been found that many industrial processes in which 
transportation lags are common can be described by a system of 
differential-difference equations, The chemical process industry offers 
many examples, 

This chapter analyzes the special case of a system subject to 
one delay, and a technique suitable for digital computation is described. 


The derivation follows a criterion developed by Koepcke (9). 
3-1 Digital solution of the matrix differential-difference equation 


Consider a dynamic system which is governed by the following 


differential-difference equation 


X(t) = AX(t) +B X(t - T) + DiU(t) + Do U(t - 7) (3.1) 
where 

X(t) = (n x 1) vector, referred to as the state vector; 

U(t) = (rf x 1) input vector, assumed constant between samples; 


i.e. U(t) = u(t) for t,<t<t 


ee KL? 


A, B = (n x n) constant coefficient matrices; and 
Di» Dy = (n x r) constant driving matrices 


Let us consider first the homogeneous part of eq. (3.1); that is 


X(t) =A X(t) + B X(t - T) (3.2) 


Taking the laplace transform of eq. (3.2), 
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SK(S) - X(0) = (A+ Be) x(s) (3.3) 
or 
x(S) = [SI - (A+ Be") x(0) (344) 
Defining =e then 
X(S) = [SE - (A+B z)] + x(0) , (3.5) 
or 
x(s) = 2 [1 - (a + Bz)/s}™* X(0) . (3.6) 
Let W= [I - Rg}, where R= A = BZ . then 
Vetiter oe en eo, (3.7) 


Therefore, one should choose an "S" large enough to ensure that 
eq. (3.7) is valid. 


Thus 


2 3 
X(S) * 2 [1 + A + BZ + ( a + (A ae + 
s s 


4 
+ A+B.) x00) (3.8) 
Ss 
Recall the facta that 
(A + Bz)? = a? + A(Bz) + (BZ)A + (Bz) 


3 


(A + Bz)? = A> + a2(Bz) + A(BZ)A + A(BZS + (BZ)AZ + 


+ (BZ)A(BZ) + (Bz)2A + (Bz)? 
(a + Bz)* = AY + a2(pz) + a2(Bz)A + A2(Bz)7+ A(BZ)A2 + 
+ A(BZ)A(BZ) + A(Bz)2A + A(BZ)> + (BZ)A> + 


+ (BZ)A7(BZ) + (BZ)A(BZ)A + (BZ)A(BZ)2 + 


+ (Bz)2a2 + (Bz)2a(Bz) + (BZ)2A + (Bz)* 
ete, 


Then, arranging by terms of equal delay, 
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2 3 44 
I 
xs) = G+4,+5+4+44——1 x0) + 
g? sg? 5” 8 
BZ A(BZ) + (BZ)A a* (Bz + ACBZ)A + BZ) Az 
+ (32 4 AG + GDA , + 
3 3 sf 


3 2 2 3 
+ AD + ADA + ABD AT + BDA dy xg) + 
5 —_ 
s 
BZ A(Bz)” +(82)A(BZ) + (3z)7a , a2(Bz)? + A(BZ)A(BZ 
a es 
s s s 
2 2 22 


+ A(BZ)"A + (BZA (BZ) + (BZ)A(BZ)A + (BZ) "A + 


3° 


+-- } X(0) + 
3 3 2 2 3 
+ Pep + A(BZ) + (BZ) A(BZ) = (B2) A(BZ) + (BZ) A + 
Ss S 


+-- ] X(0) + 


4 
oy 1 x(0) + 
Ss 
+ aeeene (3.9) 
Now, because 
ZX(t) = X(t - T) (Zee) , (3.10) 


We have 


z¥(0) = X(-T) , 
2 
z°X(0) = X(-21) , 


2>x(0) = X(-3T) , and so forth, 


Therefore, X(S) can be arranged in the following way. 


X(S) = 0(S)X(O) + 01(S)X(-T) + $2(S)X(-2T) + $3(S)X(-3T) + 


+ $,(S)X(-4T) + --- (3.11) 


where 


2 3 4 
$o(S) 2+S +h eke ------- === 

s s* gs? s* 5 

B AB + BA AB + ABA + BAZ 
#, (8) ae 

S S $ 

a?p + A“BA + ABA’ + BA? 
+ 2 EY WU 
5 
s 

p> ABS + BAB + B°A , A“B*+ABAB + ABCA + BA“B 
$2.(S) =— + 4 + 5 

S 5 $ 

BABA + B-A- 
ey ; eee nee Renee 
s 
3 3 2, 22 3 
o5(5) = Be + ABLE BAB + BOAR + BOA 
5 
S s 
4 

04 (8) © be op mee nnnnnnne nee 

s 


Rearranging terms, it follows that 


2 3 4 
A 
+ 5 } emwmemmanmen 


s s s Ss 


2 
B_, AB+ BA , ACAB + BA) + BA 
ge 4 


$9 (S) 


5 
om fH 
+ 
wl 
+ 
wl 
+ 


%, (S) 
3 Ss 


2 3 
+ AIAGB + BA) + BAT] + BAY 


5° 


2 2 2 

B + + A +B + B. 

o2(s) = + Bas Ba) , ALAB a Al, 
s s s 


B[A(AB + BA) + BA“] 
+ ; ee 


s 


3 3 2 
0,(s) = Be + AB + BAB = B(AB + BAD) | __ 


s s 


i 
| 
+ 


$,(S) 


+ 
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(3.12) 


(3.13) 


(3.14) 


(3.15) 


(3.16) 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


Let us try to find a relationship among the coefficients. With 


ie A 


Gs 
= 
= 


Figure 3.1 Array of the elements of the laplace-transformed 


transition matrices 


IT°E enSzy uz umoys Aerize vy. wWIoJ TTBYS am SpuyM Ut wept STYI 


8T 
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It is seen that the correlation among the elements (call any 


element by C ) is 


i,j 


where the arrows indicate the inmediate dependance; i.e., C,, depends on 


12 


Co and Cia etc, 


From a careful study of the array in fig. 3.1, it is found that 


Cc Ac, , BC 
ial . — batch , uv dedainl (3.22) 
gd git gd th 


where “i" is the subindex denoting row and "j" is the subindex denoting 
column. 
The following conditions should be added, in order to initialize 


a computational procedure 


Cc = 0 >0 3.23 
-1, 32: (3.23) 
Coal (3.24) 
Ce a = 0 i1>0 (3.25) 


The inverse laplace transform of eq. (3.22) yields (note: 


Le" Ait? = ayst*t 5 
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pot a % 
[€, 4] G-pt 7 (AC, 53) ST + (BC, y-1! at (3.26) 


Therefore 


c, , = (ac, , .] 3 (y-1)t + [BC Ld Gept ,(3.27) 
1,3 i ha ae a ie sa 


or 


; . [Ac, salt + (BC, 5 1-1!" 
i,j 3 : 


Changing j for j+l, eq. (3.28) takes the final form 


(3.28) 


[at]c, , + [Br]C 
- if inla4 
Cy 41 aT (3,29) 


Actually eq. (3.29) gives all coefficients without any need to 


multiply them by j 


x= 


jl ° 
This is because T has been associated with matrix A and B, and in 
order to compute any Cy yen? the initial conditions given by eqs. (3.23), 
(3.24) and (3.25) have to be considered. 


The computation of the C is done in a recursive way, as 


1,jt+l 
given by eq. (3.29). Once they are computed, they may be substituted in 
the inverse laplace transformation of eqs. (3.17), (3.18), etc., so that 
$o(T), 01(T), O2(T), .+. can be generated. The last set of matrices are 
called "plant transition matrices". 


Returning to eq. (3.11), if Pi is multiplied into both sides, 


then 
e®Sx(s) = 5(S)e°5x(0) + $(S)e°X(-T) + 02(S)e°5x(-27) + 


+ $3(s)e°5X(-3T) + —=- , (3.30) 


or 
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e*°x(s) = Oo (S)X(t) + 01(S)X(t - T) + O2(S)X(t - 27) + 
+ $3(S)X(t - 3T) +—-—. (3.31) 


Taking the inverse laplace transform of eq. (3.31), it turns out 


to be 
X(t + 1) = do(t)X(t) + 1(1) X(t - T) + O2(t)X(t - 2T) + 
+ O3(T)X(t - 3T) + -----, (3.32) 
or 
X(t +1) =) o (t)X(t - 47) (3.33) 
i= 


This is the sampled version of the homogeneous part of the 
differential-difference equation. 

Now, let us consider the addition of an input vector or forcing 
signal. 

In chapter 2, section 2=2, it was found that the digital version 


of the time-invariant matrix differential equation adopted the form 


X(K+I T) = 6(T) X(KT) + A(T) U(KT) (3.34) 
where 
o(t) = eA? (3.35) 
and 
a(t) = (eAT - 1) at (3.36) 
Although it was not demonstrated, it can be shown that 
ACT) e ane ™D, (3.37) 
j=0 
or 
A(T) = 5 ap? 1 TD. (3.38) 


j=0 
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If the terms and TD were absent, the series would be the well 


a2, 

jtl 

known matrix exponential, whose terms can be computed in a recursive way by 
- SD 

S505 441 (3.39) 


Therefore, eq. (3.38) is 
° I 
A(t) = } Coy pa? (3.40) 
j=0 
By following the same line of reasoning, the control transition 


matrices in the case of the complete differential~difference equation can 


be written as 


= tT = T 
4,1) = 2s ja + yon Fat 2 (3.41) 


and the complete difference equation is 


X(t +t) = ) [0,(t) X(t - 47) + A, (t) Ult - 47)] (3.42) 
i=0 


In resume, the digital version of 


R(t) = A X(t) + B X(t - T) + pd U(t) + Dut - T) 


is 
X(t +1) = Y [,(t) X(t = int) + A, (t) UCt - iNt)] 
i=0 
where 
T 
N= T 
fo) = 
i probes 
. 7 [at] Cc, , + (Bt] C,_, , 
L,4+1 jt+l 
Cc =I 
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CHAPTER 4 


ALGORITHMS FOR DIGITAL COMPUTATION 


This chapter presents flowcharts for the algorithms of chapters 
2 and 3, from which the computer programs were derived, They accept as 
input the coefficient matrices, the driving matrices, the initial state 
vector, and deterministic forcing vectors, As output, the computer will 
produce the state vector at the current sampling time and the set of 
transition matrices, if desired. 

Because these routines will eventually become part of Project 
ENPORT, they were designed to be used on the time-sharing system. However, 
they may be operated in the BATCH procedure without any difficulty, by 
modifying the input/output statements. 

The programs were written in the MAD language, and are listed in 


Appendix A. 
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4-1-1 TRANS 


Purpose: to compute the time response of linear time-invariant 


systems, 


Inputs: order of system (M = ); sampling time (T = ); final 
time (IF = ); number of input signals (R = ); the 


augmented A matrix and the initial state (x(1) = ). 


Outputs: the transition matrix; the current time; and the state 


of the systen. 


Remarks: main program, Subroutines called by TRANS: EXPMAT, and 


DISTUR, 
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Read from the console 


order of system (M = ) 
sampling time (T = ) 
final time (TF = ) 


Is there any disturbing 


signal ? 


Read from the console 
number of input signals 


(R=) 


Read from the console 
the augmented A matrix 

A(1,1) = --, A(2,1) = -- 
A(M,1) = -- 


Read from the console 
initial state X(1) = -- 


Save initial state 
XI(1) = X(1) 
XI(M-R) = X(M~R) 


TA is running time 


TRANS. Page 1 of 3 pages. 


Execute distur. (TA) 


XI(M-R+1) = X(M-R+1) 


XI(M) = X(M) 


Initial error 
E(1) = 0. 


E(M) = 0. 


E(I) is error vector 


TZ is time increment 


Get transition matrix 
Execute expmat. (T) 


Output to console 


Transition matrix 
EM(I,J) 


Save transition matrix 
EMP(I,J) = EM(I,J) 


Initialization 
PE(I) = 0 
Y(I) = 0 


PE is new E 


Y is new X 


TRANS. Page 2 of 3 pages 
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Compute new state 
Y(I) = Y(I) + EMP(I,J)*X(J) 
PE(I) = PE(I) + (EMP(I,J) + RIJ)*E(J) RIJ is upper bound in 
+ RIJ*X(J) remainder terms 
NorM = J |PE(I) | 
I 


Is NORM > 10°’ 2 


no 


Execute expmat. (T) 


PE(I) = 0, 
Y(I) = 0. 


Y¥(I) = Y¥(L) + EM(I,J)*XI(J) 
PE(I) = PE(I) + RIJ*XI(J) 


X(I) = Y(I) 


E(I) = PE(I) 


Output to the console 
TA, X(1) ... X(M=R) 


‘Ig TA < TF 7/29 TA = TAtTF Get new disturb. vector 


Execute distur. (TA) 
yes 
End of program (Y) 


TRANS. Page 3 of 3 pages. 


APT 


B(I,J) = A(I,J) 


B(L,J) = B(L,J)*T 


B 


a= 2° > max|B(1,J)| 
8 is the smallest positive 


integer 


B(I,J) = B(I,J)/a 


NORM = 4B = min{max(Z/B(1,J)|J, max{Z|B(I,J)|]} 
I J JI 


Initial values of K 
K = max(2*NORM, 2] 


Increment 
IN = K/2 


K is the number of 


terms of the series 


EXPMAT. Page 1 of 2 pages. 


Is RIJ < 107? | EM(1,) | 2 


LL = LL +1 


yes 
End of function 


EXPMAT. Page 2 of 2 pages. 
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3: 


qi- 
ye date 
ais 


SEM Cot 


suoreurine 


be 


furcing Sipnal vector 


called iv TRANS, 


ih? fre 
“hes 
+ Blea! 


at 


the current 
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retum 


LEO 
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4-2-1 TIMDEL, 


Purpose: to compute the time response of linear systems with 


lumped parameters and time delays. 


Inputs: order of system (M = ); sampling time (T = ); time 
delay (TD = ); final time (TF = ); number of input 
signals (R=); the A matrix; the B matrix; the D, 


matrix; the D, matrix; the initial state (X(1,1) = ). 


2 


Outputs: the plant transition matrices, the control transition 
matrices if desired; the current time; and the state 


of the systen. 


Remarks: main program. Subroutines called by TIMDEL: DELFOR, 


and PERTUR, 


Read from the console 


order of system (M = ) 
sampling time (T = ) 
time delay (TD = ) 
final time (TF = ) 


Is there any disturbing 
signal? 


Read from the console 


number of input signals 
(R=) 


REL = TD/T REL > 1, integer 


Read from the console 
the A matrix. A(1,1) = -- 
A(2,1) = == 


Read from the console 
the B matrix. B(1,1) = -- 
B(2,1) = -- 


Read from the console 
the Dy matrix. D, (1,1) =a 
D, (2,1) = 


TIMDEL. Page 1 of 3 pages. 
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Read from the console 


the Dy matrix, D, (1,1) = me 


D,(2,1) = = 


Compute set of transition 


matrices 
Execute delfor. (T) 


iio Do you wish to have the 
transition matrices ? 


Output to the console 


W is the number of 


plant transition matrices transition matrices 


EM(L), control transition computed 
matrices DELF(L) 

Read from the console 

initial state X(1,1) = --- X is (M x W) 


WI 
X(TA+T) = Y  EM(I+1) X(TA = I*REL#T) 
I=0 


= te R #0 2) 


() 


Execute pertur. (TA) 


X(TA+T) = X(TA+T) + J] DELF(I+1) U(TA - 1*REL*T) 
I=0 


TA = TA+T 


Output to the console 


TA, X(TA+T) (8) 
TIMDEL. Page 2 of 3 pages. 


he 


as 


D)(1,J) = Dy (I,J) 4T 


D,(I,J) = D,(1,J)*T 


NORM = |] Aff = Min max[(r]A(I,J) |], max[Z|A(I,J) |] 
I J Jt 


Initial value of K K is the number of terms 
K = max[{2*NORM, 2] of the series er 


Increment 
IN = K/2 
C(-1,J) = 0, J>0 


DELFOR. Page 1 of 2 pages. 
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c(1,s+1) = CALCD + (B1¢d-1,y) 


J+1 


K 
EM(I+1) = J) C(1,J) 
J=0 


= C(I gJ) *TAD 
G(L,I) = }) Sete L 


J=I 


yes 


Is I #0? 


n 


8 
Qo 


RIJ is upper bound in 


remainder terms of Pe 


Ie RIJ < 10°’ |EM(x)| 2 


yes 


Ie |lemc1+1)|| < 10772 


yes 


Function return 


End of function 


DELFOR. Page 2 of 2 pages. 


Aemarks | 


Bue 


PERTUR 


ompute the forcing Ssignui vector at tne current 


. lhe program has to Keep track of the past. 


gutine called by ITMbDEL. 


41 


¥ 


waa! 


Hae 


CHAPTER 5 


SOLUTION TO SAMPLE PROBLEMS 


This chapter describes a set of sample problems which were 
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selected because they represent typical applications of the two simulators. 


They are intended to show the use of the state variable diagram, and also 


to show the accuracy of the methods, 


Dd-1 Test problem for the simulation of dynamic systems without delay 


Example 5-1, Although this example may represent a great number 


of physical processes, it was selected purely from the mathematical point 


of view, The same problem was run by Liou (11). 


Given 


and 


Obtain X(nT) using T = 0.1 Min. 


The reported solution by Liou and the one obtained by the 


simulator are 


(5.1) 


(5.2) 


loadgo trans expmat distur 


W 1010.4 
EXECUTION. 

GIVE ORDER OF SY 
SAMPLING TIME (T 
m=3,t=.1,tf=2.* 


STEM (M 
), 


) 


FINAL TIME (TF = ) 


IS THERE ANY DISTURBING SIGNAL 


no 


GIVE THE A MATRIX (AC1,1)=--,A(2,1)=--) 


a(1,1)=20.,1.,0.* 
a(2,1)=0.,0.,1.% 
a(3,1)=-.75,-2.7 


GIVE 
x(1)22.,°2.5,3.7 


5,73.* 


Se 


INITIAL STATE (X(1)=--) 


TERMS OF THE MATRIX EXPONENTIAL 


EM 1 
EM( 1 
EM ( 1 
EM( 2 
EM( 2 
EM( 2 
EM( 3 
EM( 3 
EM( 


TIME 


-10 
20 
- 30 
40 
50 
60 
-70 
80 
90 
1.00 
1.10 
1.20 
1.30 
1.40 
1.50 
1.60 
1.70 
1.80 
1.90 
2.00 
2.10 
END OF EXECUTION 
TO CONTINUE, 


, 1) 
’ 2) 
, 3) 
, 1) 
’ 2) 
, 3) 
’ 1) 
? 2) 
’ 3) 


X(1) 


-176781E 
«L56774E 
~139515E 
-124603E 
-111700E 
-LOQS515E 
-907978E 
»823379E 
»749538E 
-68K4911E 
-628178E 
-578215E 
-534062E 
-494901E 
-4600354E 
«428868E 
-400894E 
-375681E 
-352861E 
-332118E 
-313185E 


AND PRINT AN ASTERISK 


-999884E 00 
-995717E-01 
~452513E-02 
-.339385E-02 
-987440E 00 
- 859963E-01 
-.644972E-01 
-.23988KE 00 
»729451E 00 


X(2) 
ol -.215290E 
ol -.185614E 
ol -.160242E 
01 -.138548E 
01 -.119997E 
ol -.104131E 
00 -.905571E 
00 -.7894K13E 
00 ~.689964E 
00 -.604775E 
00 -.5317535E 
00 -.469107E 
00 -.415312E 
00 -.369064E 
00 -.329250E 
00 -.294921E 
00 -.265268E 
00 -.239602E 
00 -.217334E 
00 -.197965E 
00 -.181068E 


GO TO THE TOP OF A NEW PAGE 


ee er as erence en een an eer ane eee 


01 
01 
o1 
01 
ol 
01 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 


V(t) = A (te) 


X(1) 
V(t) =/X(2) 
X(3) 


0 1 


A=; Q 0) 
-.75 -2.75 


2 
WO) = |-2.5 
3.75) 


X(3) 


.320616E 
~274116E 
-234368E 
- 200401E 
2171382E 
» LK6596E 
© 125431E 
~107362E 
~9I9KIBE 
-787838E 
-675590E 
-579853E 
~498212E 
-428602E 
~369257E 
-318666E 
~275537E 
-238768E 
- 207415E 
-180676E 
»157862E 


Figure 5.1 Console transaction for example 5.1 


ol 
01 
ol 
01 
ol 
Ql 
ol 
01 
00 
00 
00 
00 
00 
00 
00 
00 
00 
ao 
00 
00 
00 
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X(1), X(2), X(3) 


X,(0) X,(0) x, (0) 


1.0 1.1 1.2 1.3 1.41.5 1.6 1,7 1.8 1.9 210-251 


Po oS t (min) 


7 Figure 5.2 Response curves of system for example 5.1 


94 
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A liquid stream enters tank 1 (figure 5.3) at a volumetric flow 
rate F cfm and contains reactant A at a concentration of Cy moles Aset>, 
Reactant A decomposes in the tanks according to the irreversible chemical 
reaction, 

Aa-—— B 
The reaction is first order and proceeds at a rate 
r@ke 
where 
r= moles A decomposing/ (ft) (time) 
c ™ concentration of A, moles A/tt? 
k = velocity constant, a function of temperature 

The reaction is to be carried out in a series of two stirred 
tanks, The tanks are maintained at different temperatures. The temperature 
in tank 2 is to be greater than the temperature in tank 1, with the result 
that kos the velocity constant in tank 2, is greater than in tank 1, k,. 
Changes in physical properties due to chemical reaction are neglected, 

The purpose of the control system is to maintain Cos the 
concentration of A leaving tank 2, at some desired value in spite of 
variation in inlet concentration Coe This will be accomplished by adding a 
stream of pure A to tank 1 through a control valve. 

Further assumptions are that the control valve and the measuring 
element have no dynamics, and that the controller exert proportional action 
on the process, 

A portion of the liquid leaving tank 2 is continuously withdrawn 
through a sample line. The measuring element is remotely located from the 
process, because rigid ambient conditions must be maintained for accurate 


concentration measurements. The sample line can be represented by a 
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transportation lag. 


ontroller 


composition 


heating coils 
Figure 5.3 


Control of a stirred-tank chemical reactor 


The following data is assumed to apply to the system 
Molecular weight of A = 100 1lb/lb mole 
7 0.8 lb mole/ft” 
Cj, = O61 1b mole A/ft® 
F = 100 cfm 
m = 1.0 lb mole/min 


k, © 1/6 ain? 


k, = 2/3 min? 


v = 300 ft? 


Valve sensitivity k= 1/6 cfm/psi 


Measuring device sensitivity 


ke = 100 in. pen travel/(1b mole/ft”) 


Time delay in sample line = T 


The overall block diagram which the authors propose is 
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Figure 5.4 


Block diagram for a chemical reactor 


control system 


It is assumed that the inlet concentration fy does not change 
with time. 
As was discussed in chapter 2, the state variable diagram can 


be obtained in three ways. Direct programming will be used in this case. 


With this purpose, let it be called cy the input to the lag term ir 
and c its output in figure 5.4, then 
c 1 
es 62) 
B 
or 
-1 
c Ss 
aoe 65 Tesi: (5.4) 
B 
Eq. (5.4) can be written as 
-1 
ec=,.5S8 B, (5.5) 
where 
c. 
E ad B ~1 (5.6) 
1+,5S 
Transposing 
1 


Esc, 7 65 s" gE, (5.7) 
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By a similar procedure 


= = oy , (5.8) 
or 
cB g7l 
pa Sore 5 (5.9) 
S,  yag7} 
-1 
SS 2 § E (5.10) 
where 
€ 
z= —; (5.11) 
148 
Transposing 
a1 
Rt eg ee (5.12) 


The state variable diagram follows from eqs. (5.5), (5.7) and 


eqs. (5.10), (5.12), and is shown in figure 5.5. 


Figure 5.5 
State variable diagram for a chemical reactor 


control systen 


The notation in figure 5.5 has been changed slightly. This is in 
order to follow the same symbolism given in the previous chapters. 


In figure 5.5 the state variables ara X, and Xe The differential- 
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difference equations for the atate variables are readly obtained by 


inspection of the diagram, That is, 


Xp 7.5 X) + 05 xX, (5.13) 
X,=-X, + K U-K X(t - T) (5.14) 
Therefore the matrix differential-difference equation is 
e -.5 ° Q 0 0 
X(t) = X(t) + X(t - T) + U(t) (5.15) 
Q -1 -K (¢) K 
where 
X(t) 
X(t) = 
X, (t) 


From this equation, it is seen that the coefficient matrices and 


driving matrices are 


-.5 45 
A= (5.16) 
0 -1 
0 ft) 
B= (5.17) 
-K 0 
0 
D Cd (5.18) 
z K 
0 
D, = (5.19) 
2 0 


Five numerical examples were run using this system. These are 


summarized as follows. 


Example 5.2.1. Overall forward gain K= 5,24. We assume a time delay 
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equal to zero. A unit step is the input and all initial conditions are 


zero, The matrix differential equation is 


-.5 oe) ce) 


ve) = [75.24 -1 5.24) V(t) (5.20) 
) ) 0 
where 
X(t) 
M(t) = |X, (t) (5.21) 
UCt) 


Example 5.2.2. Overall forward gain K = 5,24, and time delay = 
»5> Mins. Same conditions of the state were taken, The matrix differential- 


difference equation is 


. -.5 65 0) 0 ) 
X(t) = X(t) + X(t - .5) + U(t) + 
) -1 “5.24 0 5424 


) 
+ fe - .5) (5.22) 
0 


Example 5,2.3. Overall forward gain K = 1.85. Time delay is zero. 


The remaining conditions are the same. The state equation is 


-=5 5 O 
W(t) = {-1.85  -1 1.85 v(t) (5.23) 


0 Q 0 


Example 5.2.4. Overall forward gain K = 1.85. Time delay = .5 min, 


Unit step and zero initial conditions are assumed. The state equation is 
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7 -.5 25 i) 0) ) 
X(t) = X(t) + X(t - .5) + u(t) + 
0 -1 -1.85 0 1,85 
) 
+ U(t - .5) (5,24) 
0 


Example 5.2.5. This is the same as example 5.2.4, with the 
exception of the time delay, which is taken equal to 1 min.. The state 


equation is 


. -.5 5 i?) 0 0 
X(t) = X(t) + X(t - 1) + 3) 
0 -1 ~1.85 0 1.85 


0 
+ | u(t = 1) (5,25) 
0 


All five examples with the input/output information and the 
response curves, are shown in figures 5.6 to 5.15. 

The interested reader should compare the responses of the three 
cases with delay with those given by Coughanowr and Koppel on page 467 of 


reference (4). 
5-3 Test problem for the simulation of dynamic systems with delays 


The eighth example was run in order to check the accuracy of 
evaluation of the set of transition matrices. This example is discussed by 
Koepcke (9). 

The problem is described as an unstable process which is governed 


by 
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* 


GIVE ORDER OF SYSTEM (M 


SAMPLING TIME (T = ), FINAL TIME (TF = ) 


in=2,t=.5,tfell.* 


IS THERE ANY DISTURBING SIGNAL 


yes 


GIVE NUMBER OF 


r=l* 


GIVE THE A MATRIX (A(1,1) 


a(1,1)=-.5,.5,0.* 
a(2,1)5-5.24,-1.,5.24* 
a(3,1)=0.,0.,0.* 


GIVE INITIAL STATE (X(1)=--) 


x(1)=0.,0.* 


) 


INPUT SIGNALS (R = ) 


=--,A(2,1)=--) 


TERMS OF THE MATRIX EXPONENTIAL 


EM( 
EM( 
EM( 
EM( 
EM( 
EM( 
EM ( 
EM ( 
EM( 


TIME 


50 
1.00 
1.50 
2.00 
2.50 
3.00 
3.50 
4,00 
4.50 
5.00 
5.50 
6.00 
6.50 
7.00 
7.50 
8.00 
8.50 
9.00 
39.50 

10.00 
10.50 
11.900 


1, 1) 
1, 2) 
1, 3) 
2, 1) 
2, 2) 
2, 3) 
3, 1) 
3 2) 
3, 3) 


une btw uw wonwu 


X(1) 


~243387E 
»0605063E 
~354087E 
-103181E 
~969739E 
»873563E 
~810740E 
»795981E 
~811516E 
-833372E 
-846973E 
~849679E 
»845848E 
~840898E 
-837966E 
»837495E 
~838429E 
»839546E 
«840175E 
~840250E 
»840025E 
~859774E 


END OF EXECUTION 
GO TO THE TOP OF A NEW PAGE 
AND PRINT AN ASTERISK 


TU CONTINUE, 


Figure 5.6 Console transaction 


-5560765E 00 
-154089E 00 
-243387E 00 
-161485E 01 
-401987E 00 
~185824E O1 
-900000E 00 
-O00000E 00 


1.000000E 00 


00 
00 
00 
01 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 


X(2) 


~185824E 
~221219E 
~167353E 
-990269E 
~590102E 
~529468E 
-660403E 
~814488E 
-900262E 
~909654E 
~878156E 
~843503E 
~825210E 
~824044E 
-831568E 
~839327E 
~843206E 
~843258E 
~S41475E 
- 8397K3E 
-838925E 
- 833960E 


V(t) = A V(t) 


X(1) 
M(t) = 1X(2) 
u(1) 


7.5 
A= 1[-5.24 
0 


for example 5.2.1 


) 
-1 
0 


0 
5.2 
0 


Figure 5.7 Response curves of system for example 5.2.1 
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loadgo timdel delfor pertur 
Ww 1039.4 

EXECUTION. 

GIVE ORDER OF SYSTEM (M = ) 


i} 
DESIRED SAMPLING TIME (T = ) i 
TIME DELAY (TD = ), FINAL TIME (TF = ) | 
m=2,t=.5,td=.5,tf=15.* 1 
', 
1S THERE ANY DISTURBING SIGNAL bake) = A X(t) + B X(t = 5) 
es 
H + D,U(t) + DUE ~ .5) 
NUMBER OF INPUT SIGNALS (R = ) i 
r=l* H 
i 
GIVE THE A MATRIX (A(1,1)=--,A(2,1)=--) H 
a(1,1)=-.5,.5* Tee ee “A 
a(2,1)=0.,-1.* 1 0 “1 
t 
t 
GIVE THE B MATRIX (B(1,1)=--,8(2,1)=--) ! 
b(1,1)=0.,0.* ate a) ) 
b(2,1)=-5.24,0.* ! -5.24 0 
t 
t 
GIVE THE Dl MATRIX (D1(1,1)=--,01(2,1)=--) | 
dl(1,1)=0.* fe cee 2 | 
d1(2,1)=5,24* aes 5424 
t 
i 
GIVE THE D2 MATRIX (02(¢1,1)=--,02(2,1)=--) } _ 
d2(1,1)=0.* ‘— 0 
d2(2,1)=0.* 1 "2 0 
DO YOU WISH TO HAVE THE TRANSITION MATRICES 
yes 
TRANSFER MATRIX PHI ( 0) TRANSFER MATRIX DELTA 0) 
EM ( 1, 1) =  .778801£ 00 DEL( ly 1) = .256388E 00 
EM( 1; 2) =  .172270E 00 
EM( 2, 1) = .000000E 00 DEL( 2, 1) =  .206178E 01 
EM( 2, 2) = .606531EF 00 
TRANSFER MATRIX PHI ( 1) TRANSFER MATRIX DELTA( 1) 
EM( 1, 1) = -.235067E 00 DEL( 1, 1) = -.132818E-01 
EM( 1, 2) = =.187866E-01 
EM( 2:5 1) = -.180539E 01 DEL(— 2, 1) = -.210165E 00 
EM( 2, 2) = -.216281E 00 
TRANSFER MATRIX  PHI( 2) TRANSFER MATRIX DELTA(C 2) 
EM( 1, 1) = .126127£-01 DEL ly 1) = = .283552E-03 
EM( 1, 2) 5 .614986E-03 
EM( 2, 1) = .196883E 00 DEL( 2, 1) =  .672861E-02 
EM( 25 2) = .119977E-01 
TRANSFER MATRIX  PHI( 3) TRANSFER MATRIX DELTA( 3) 
EM( 1, 1) = -.273338E-03 DEL( 1; 1) = -.327556E-05 
EM( 1, 2) = =-.958848E-05 
EM( 2, 1) = -.64450G6E-02 DEL( 2, 1) = -.103763E-05 
EM( 2, 2) = -.263750E-03 
TRANSFER MATRIX  PHI( 4) TRANSFER MATRIX DELTAC 4) 
EM( 1, 1) =  .318384E-05 DELC 1, 1) =  .236515E-07 
EM( 1, 2) = .872148E-07 
EM ( 2, 1) = .100487E-03 DELC 25 1) = .937662E-06 
EM ( 2, 2) =  .309662E-05 
TRANSFER MATRIX PHI( 5) TRANSFER MATRIX DELTAC 5) 
EM( 1, 1) = ~.231099€-07 DEL( 1; 1) = -.116723£-03 
EM( 1, 2) = ~.5192686-09 
EM( 2, 1) = -.914011E-06 DELC 23 1) = -.555865E-08 


EM( 2, 2) -.225906E-07 
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TRANSFER MATRIX  PHI( 6) TRANSFER MATRIX DELTA(C 6) 
Em ( 1, 1) = »114463E-09 DEL 1, 1) = -418403E-12 
EM( 1, 2) = -218008E-11 
EM ( 2, 1) = »544192E-08 DEL( 2, 1) = »232657E-10 
EM( 2, 2) = ©112283E-09 
GIVE THE INITIAL STATE (X(1,1)=---) 
x(1,1)=0.,0.* 
TIME X(1) X(2) 

50 -2564E 00 -2062E O01 
1.00 ~7980E 00 -3102E O1 
1.50 -1300E 01 ~2831E 01 
2.00 -1502E 01 -1539E O01 
2.50 -1332E 01 -2406E-01 
3.00 -9204E 00 -.8884E 00 
3.50 -5132E 00 -.7847E 00 
4.00 -3246E 00 -1653F 00 
4.50 -4295E 00 -1364E O01 
5.00 -7391E 00 -2150E O01 
5.50 »1067E O1 ~2155E O01 
6.00 -1237E O01 -1465E 01 
6.50 -LI79E O1 ~5227E 00 
7.00 ~9475E 00 -.1450E 00 
7.50 -6857E 00 -.2169E 00 
8.00 »S349E 00 -2774E 00 
8.50 -5622E 00 -1013E O1 
9.00 -7330E 00 »1573E O12 
9,50 -9408E 00 ~1683E 01 

10.00 »1072F O01 »1335E O1 
10.50 »1065E Ol -7646E 00 
11.00 -SKOLE 00 ~2988E 00 
11.50 -7766E 00 -1715E 00 
12.00 -6646E 00 «411I5E 090 
12.50 -6579E 00 -8506E 00 
13.00 «7480E 00 -1234E 01 
13.50 -8763E 00 -1366E O01 
14.00 -9708E 00 ~1205EF 01 
14.50 -9854E 00 -8692E 00 
15.00 ~9214E 00 -5560E 00 


END OF EXECUTION 
TO CONTINUE, GO TO THE TOP OF A NEW PAGE 
AND PRINT AN ASTERISK 


Figure 5,8 Console transaction for example 5.2.2 


Us unit step 


Figure 5.9 Response curves of system for example 5.2.2 


BS 


* 


GIVE ORDER OF SYSTEM (M 


SAMPLING TIME (T = ), 


in=2,t=.5,tfell.* 


IS THERE ANY DISTURBING SIGNAL 


yes 


GIVE NUMBER OF 


r=l* 


GIVE THE A MATRIX (A(1,1)=--,A(2,1)=--) 


a(1,1)=#-.5,.5,0.* 
a(2,1)2-1.85,-1.,1.85* 
a(3,1)20.,0.,0.* 


GIVE INITIAL STATE (X(1)=--) 


x(1)=0.,0.* 


) 


FINAL TIME (TF = ) 


INPUT SIGNALS (R = ) 


TERMS OF THE MATRIX EXPONENTIAL 


EM( 
EM( 
EM( 
EM( 
EM( 
EM( 
EM( 
EM( 
EM( 


TIME 


20 
1.00 
1.50 
2.00 
2.50 
3.00 
3.50 
4.00 
&.50 
5.00 
5.50 
6.00 
6.50 
7.00 
7.50 
8.00 
8.50 
9.00 
9.50 

10.00 
10.50 
11.00 


1, 1) 
1, 2) 
1, 3) 
2, 1) 
2, 2) 
2, 3) 
3, 1) 
3, 2) 
3, 3) 


nunituaun#a 


X(1) 


-888757E- 


.267189E 
»444358E 
-577874E 
-658281E 
-694033E 
-699993E 
»-690429E 
»-675860E 
~-662472E 
~652899E 
-647458E 
-645293E 
-645202E 
-646113E 
~647276E 
-O48274KE 
-6489535E 
~G4931KE 
-649438E 
»-G49420E 
~6493359E 


END OF EXECUTION 
GO TO THE TOP OF A NEW PAGE 
AND PRINT AN ASTERISK 


TO CONTINUE, 


Figure 5.10 Console transaction for example 5.2.3 


1 


01 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 


-697370E 00 
-165714E 00 
-888757E-01 
«613141E 00 
-531656E 00 
-702016E 00 
-QO0000E 00 
-QOO000E 00 
-OO0000E 00 


X(2) 


-702016E 
~102075E 
-108088E 
-100422E 
~881598E 
-767104E 
~684312E 
-636641E 
-617160E 
-615736E 
-623187E 
-633019E 
-G41I581E 
-647461E 
-650643E 
-651776E 
-651666E 
-650995E 
-650222E 
-649590E 
-649178E 
-648970E 


v(t) = A V(t) 


u(1) 


X(1) 
v(t) =|x(2) 


-.5 
A«- ~1.8 
0 


X(0) -| 


00 
01 
Ol 
01 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 
00 


5 


0 
0 


5 i¢) 
-1 1.8 
0 0 
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Us unit step 


Figure 5.11 Response curves of system for example 5.2.3 


09 


* 

GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY (TD = ), FINAL TIME (TF = ) 
m=2,t=.5,td=.5,tf=15.* 


IS THERE ANY DISTURBING SIGNAL 
yes 


NUMBER OF INPUT SIGNALS (R = ) 
r=l* 


GIVE THE A MATRIX (A(1,1)=--,A(2,1)=--) 
a(1,1)=-.5,.5* 
a(2,1)0.,-1.* 


GIVE THE B MATRIX (B8(1,1)=--,B(2,1)=--) 
b(1,1)=0.,0.* 
b(2,1)=-1.85,0.* 


GIVE THE D1 MATRIX (D1(1,1)=--,01(2,1)=--) 
d1(1,1)50.* 
d1(2,1)=1.85* 


GIVE THE D2 MATRIX (02(1,1)=--,D2(2,1)=--) 
d2(1,1)=0.* 

d2(2,1)=0.* 

DO YOU WISH TO HAVE THE TRANSITION MATRICES 
no 


ws ce ce cee en ce ere ee ee a ees a a ee ee OT OD oO AR eS a TO 


GIVE THE INITIAL STATE (X(1,1)=---) 
x€1,1)20.,0.* 
TIME X(1) X(2) 

50 ~9052E-01 ~7279E 00 
1.00 -2848E 00 -1l143E 01 
1.50 -4952E 00 -1282E 01 
2.00 -6644E 00 -1216E 01 
2.50 -7664E 00 -1034E O1 
3.00 -8015E 00 -8266E 00 
3.50 -7862E 00 -6539E 00 
4,00 ~74H31E 00 -S448E 00 
4.50 -6932E 00 -5021E 00 
5.00 -6512E 00 -5114E 00 
5.50 -6245E 00 -5508E 00 
6.00 -6138E 00 -5995E 00 
6.50 -6158E 00 -6421E 00 
7.00 -6251E 00 -6704E 00 
7.50 -G369E 00 -6829E 00 
8.00 ~6472E 00 -6825E 00 
8.50 ~6541E 00 -6740E 00 
9.00 -6572E 00 -6627E 00 
9.50 -6572E 00 -6523E 00 

10.00 -6552E 00 -6450E 00 
10.50 -6525E 00 -6414E 00 
11.00 -6499E 00 -6411E 00 
11.50 -6482E 00 -6429E 00 
12.00 -6473E 00 -6455E 00 
12.50 -6472E 00 -G480E 00 
13.00 -6476E 00 -6499E 00 
13.50 -6482E 00 -6508E 00 
14.00 ~G488E 00 -6510E 00 
14.50 -G4I3E 00 -6507E 00 


15.00 =6495E 00 -G501E 00 
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X(t) =A X(t) +B X(t - .5) 


+ D,U(t) + D,U(t - .5) 


0 
eS las 
0 
i ioe 
0 
te [ 0 


Figure 5,12 


5.2.4 


Console transaction for example 


(4) 


x(1) 
X(2) 


2 4 6 8 10 12 14 t (min) 


Figure 5.13 Response curves of system for example 5.2.4 
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loadgo timdel delfor pertur 

W 1250.1 

EXECUTION. 

GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY (TD = ), FINAL THAME (TF = ) 
in=2,t=.5,td=1.,tf=15.* 


IS THERE ANY DISTURBING SIGNAL X(t) = A X(t) +B X(t - 1) 
yes = ve 


NUMBER OF INPUT SIGNALS (K = ) PPiNCe) “FD UCe =) 


r=l* 


i] 

1 

i 

1 

{ 

1 

1 

i 

1 

i 

i 

1 

1 

{ 

H 

GIVE THE A MATRIX (AC1,1)=--,A(2,1)2--) 1 
a(l,1)=-.5,.5* 3 Ss 
a(2,1)20.,-1.* Af [3 S| 

GIVE THE B MATRIX (B(1,1)=--,8(2,1)=--) H 

H 

i 

J 

1 

1 

1 

i} 

1 

! 

i 

1 

1 

i] 

1 

i] 

1 

t 

BS) 


b(1,1)=0.,0.* 
b(2,1)=-1.85,0.* 


GIVE THE D1 MATRIX (D1(1,1)=--,01(2,1)=--) 
di(1,1)=0.* 
d1(2,1)21.85* 


GIVE THE D2 MATRIX (D2(1,1)=--,02(2,1)=--) 
d2(1,1)=0.* 
d2(2,1)=0.* 
DO YOU WISH TO HAVE THE TRANSITION MATRICE 


yes 

TRANSFER MATRIX PHI( 0) TRANSFER MATRIX DELTA( 0) 
EM( 1, 1) = -778801E 00 DELC 1, 1) = -905188E-01 
EM( 1, 2) = ~172270E 00 

EM ( 2, 1) = -OO0000E 00 DELC 2, 1 = »727918E 00 
EM ( 2, 2) = -606531E 00 

TRANSFER MATRIX  PHI( 1) TRANSFER MATRIX DELTAC 1) 
EM( 1, 1) = -.829913E-01 DELC 1, 1) = -.165553E-02 
EM ( 1, 2) = ~.663267E-02 

EM( 2, 1) = -.637399E 00 DEL 2, 1) = -.261964E-01 
EM( 2, 2) = -.763586E-01 

TRANSFER MATRIX  PHI( 2) TRANSFER MATRIX DELTAC 2) 
EM( 1, 1) = ~157213E-02 DEL( 1, 1) = ~124783E-08 
EM( 1, 2) = . 766560E-04 

EM ( 2, 1) = » 245409E-01 DELC 2, 1) = ~296106E-03 
EM( 2, 2) = ~149548E-02 

TRANSFER MATRIX PHIC 3) TRANSFER MATRIX DELTAC 3) 
EM( 1, 1) = -.120288E-04 DELC 1, 1) = -.508918E-07 
EM( 1, 2) = -,421960E-06 

EM ( 2, 1) = -.283627E-03 DELC 2, 1) = -.161214E-05 
EM( 2, 2) = -.116068E-04 

TRANSFER MATRIX PHI( 4) TRANSFER MATRIX DELTAC 4) 
EM(¢ 1, 1) = -494666E-07 DEL( 1, 1) = ~129736E-09 
EM( 1, 2) = ~135504E-08 

EM( 2, 1) = -156125E-05 DELC 2, 1) = -514338E-08 
EM( 2, 2) 5 ~481116E-07 

TRANSFER MATRIX PHI¢ 5) TRANSFER MATRIX DELTAC 5) 
EM( 1, 1) = -.126765£-09 DEL( 1, 1) = -.226048E-12 
EM( 1, 2) = -.284835E-11 

EM( 2, 1) = -.501365E-08 DEL( 2, 1) = -.107649E-10 
EM( 2, 2) = ~-.123917E-09 
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GIvE THE INITIAL STATE (X(1,1)=---) 
x€1,1)=0.,0.* 
TIME X¢1) X(2) 
50 -9052E-01 ~7279E 
1.00 -2864E 00 ~1169E 
1.50 -S154E 00 ~141I1E 
2.00 -7194E 00 -LS44E 
2.50 -8664E 00 ~1306E 
3.00 -9369E 00 -1063E 
3.50 -9328E 00 ~ 7864E 
4.00 -8712E 00 ~5417E 
4.50 -7771E 00 »3720E 
5.00 -6770E 00 -2960E 
5.50 -5928E 00 ~3093E 
6.00 -5384E 00 ~3897E 
6.50 ~5185E 00 -5062E 
7.00 ~5299E 00 -6269E 
7.50 ~5634E 00 ~7262E 
8.00 -6073E 00 »7881E 
8.50 -6503E 00 - 8080E 
9.00 -6838E 00 ~7909E 
9.50 -7029E 00 -7482E 
10.00 -7068E 00 -6944KE 
10.50 -6980E 00 »6429E 
11.00 -6810E 00 -6038E 
11.50 -G611E 00 ~5825E 
12.00 -6430E 00 ~5795E 
12.50 -6301E 00 ~591S5E 
13.00 -6239E 00 -6128E 
13.50 -6243E 00 ~6370E 
14.00 -6296E 00 ~-6584E 
14.50 -6379E 00 ~6733E 
15.00 -6466E 00 -6800E 


END OF EXECUTION 
GO TO THE TOP OF A NEW 
AND PRINT AN ASTERISK 


TO CONTINUE, 


Figure 5.14 Console transaction for example 5.2.5 


X(1), X(2) 


+++ +--+ 
2 4 6 8 10 12 14 16 ¢t (min) 


Figure 5.15 Response curves of system for example 5.2.5 
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. ) 1 2 ) 0 
X(t) = X(t) + X(t = T) + U(t - T) (5.26) 
-1 0 0 - 1 1 


Te is assumed the sampling time equal to the time delay. That is, 
Tete —b min (5.27) 


Koepcke reported the following results of the plant transition 


matrices and the control transition matrices: 


«7071068 7071068 «00000 
% = Ao = 
=. 7071068 «7071068 00000 
«1338340 20277680 2928932 
% = 4, = 
-.0277680 -.0782980 «7071068 
20109582 20022524 20075873 
% = 42 * 
=,0022524 0026278 -.0308106 
20005903 20000742 20004532 
5 = A3 = 
=,0000742 -,0000854 20007349 
«0000236 0000026 0000119 
, s Ay a 
-.0000026 20000013 =.0000165 
0000008 20000001 -0000003 
a = 45 = 
=-,.0000001 -,0000000 «0000002 


The time response of the syatem was obtained assuming a step 


input and zero initial conditions for the integrators. 


The solution is depicted in figures 5.18 and 5.19. 


In a similar way, this same example was tested assuming no lags 


in the system, that is 


02 1 0 
wt) = fen =. af vee) (5.28) 
0 0 0 
where 
X, (¢) 
M(t) = 1x, (e) (5429) 
U(t) 


The evaluation of the state is shown in figures 5.16 and 5.17. 

It is interesting to compare the transient response in both cases. 
As it can be seen in the plots (figures 5.17 and 5,19), the case with 
delay is something less unstable than the linear one with delay equal to 


zero. 
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* 


GIVE ORDER OF SYSTEM (M = ) 


SAMPLING TIME (T = ), FINAL TIME (TF = ) | 
m=2,t=.5,tf=l5.* 

1 
1S THERE ANY DISTURBING SIGNAL V(t) = A V(t) 
yes ! 

H X(1) 
GIVE NUMBER OF INPUT SIGNALS (R = ) ¥ v(t) =| x(2) 
r=l* = U(1) 
GIVE THE A MATRIX (A(1,1)=--,A(2,1)=--) | 
a(1,1)=.2,1.,0.* i 201 
a(2,1)=-1.,-.1,1.* Vase | -1  -.1 
a(3,1)=0.,0.,0.* 0 0 

{ 
GIVE INITIAL STATE (X(1)=--) | 0 
x(1)=0.,0.* H x(0) “3 | 


TERMS OF THE MATRIX EXPONENTIAL 
EM ( 1, 1) -976370E 00 
EM( 1, 2) -492031E 00 


EM (¢ 1, 3) ~124527E 00 
EM ( 2, 1) -.492031E 00 
EM( 2, 2) -828760E 00 
EM( 2, 3) ~467126E 00 


EM( 3, 1) 
EM ( 3, 2) 
EM(¢ 3, 3) 


-O00000E 00 
-OO0000E 00 
1.000000E 00 


one bt tt we Roo 


TIME X¢1) X(2) 

50 ~124527E 00 -467126E 00 
1.00 -475952E 00 -7929390E 00 
1.50 -979407E 00 -890141E 00 
2.00 -151877E 01 »722940E 00 
2,50 -196311E O1 -318989E 00 
3.00 ~219820E 01 -.234422E 00 
3.50 »215544E O1 -.808739E 00 
4.00 -183111E O01 -.126367E O01 
4.50 -129060E Ol -.148112E 01 
5.00 -655877E 00 -.139538E O01 
5.50 -783336E-01 -.101203E O01 
6.00 -.296938E 00 ~.410143E 00 
6.50 -.367197E 00 -273318E 00 
7.00 -.995124E-01 -874313E 00 
7.50 -457555E 00 -124068E 01 
8.00 -118173E O1 -127022E O01 
8.50 -190332E 01 -938392E 00 
9.00 ~-244459E 01 -308336E 00 
9.50 -266306E O1 -.480150E 00 

10.00 »-248841E O1 -.124111E O01 
10.50 ~194347E 01 -.178583E 01 
11.00 2114338E 01 -.196915E 01 
11.50 -272011E 00 -.172740E O1 


12.00 ~.459826E 00 


-109832E O1 


“anule 


2.416 


u 


-,2)s3u"t 20 
et LZ 7226 0 
»leoglZe YQ 
J1LSa4e te 3 

? 
1 


BO G5 te 
mw 
Teo 
a 
ben be eet he 


sd 
cn 


Console transaction tor example 5.3 


when time delay = 4% 


Uz unit step 


U 


Figure 5.17 


Responge curves of system for example 5.3 when time delay = 0 


od 


rat 


* 

GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY (TD = ), FINAL TIME (TF = ) 
im=2,t=.7853982,td=.7853982,tf=20.* 


{S THERE ANY DISTURBING SIGNAL X(t) =A X(t) +B X(t - Dp 
yes — ii) a 


NUMBER OF INPUT SIGNALS (R = ) + D,U(t) + D,UCt -? 


r=l* 


GIVE THE A MATRIX (AC1,1)=--,AC2,1)=--) 
a(1,1)=0.,1.* 
a(2,1)+-1.,0.+* 


GIVE THE B MATRIX (8(¢1,1)=--,8(2,1)=--~) 


1 
1 
‘ 
1 
1 
1 
1 
1 
1 
! 
1 
1 
| 
! 
1 
' 
\ 
b(1,1)=.2,0.* | ae 2 0 
b(2,1)=0.,-.1* ! 0 1 
i 
1 
H 
1 
1 
{ 
1 
i 
1 
1 
1 
1 
1 
| 


GIVE THE D1 MATRIX (D1(1,1)=--,01(2,1)=--) 
di(1,1)=0.* 
d1(2,1)=0.* 


GIVE THE D2 MATRIX (02(1,1)=--,D2(2,1)=--) 
d2(1,1)=0.* 

d2(2,1)=1.* 

DO YOU WISH TO HAVE THE TRANSITION MATRICES 


yes 

TRANSFER MATRIX  PHI( 0) TRANSFER MATRIX DELTA(C 0) 
EM( 1, 1) = ~707107E 00 DELC 1, 1) = »OO00000E 00 
EM( 1, 2) = .707107E 00 

EM( 2, 1) = -.707107E 00 DEL ea 1) = ~O00000E 00 
EM( 2, 2) = .707107E 00 

TRANSFER MATRIX  PHI( 1) TRANSFER MATRIX DELTAC 1) 

EM ( 1, 1) = ~133834E 00 DEL( 1, 1) = ~292893E 00 
EM( Ly 2) 5 ~277680E-01 

EM ( 23 1) = -.277680E-01 DEL( 2, 1) = ~707107E 00 
EM ( 2, 2) = -,.782980E-01 

TRANSFER MATRIX PHI( 2) TRANSFER MATRIX DELTAC 2) 
EM( 1, 1) = ~109582E-01 DELC 1, 1) = ~758732E-02 
EM( 1, 2) = ~225237E-02 

EM( 2, 1) = =-.225237E-02 DEL( 2, 1) = ~.308106E-01 
EM( 2, 2) = ~262783E-02 

TRANSFER MATRIX PHI ( 3) TRANSFER MATRIX DELTA( >) 
EM( 1, 1) = -590343E-03 DEL( 1, 1) = ©453237E-03 
EM( 1, 2) = ~741765E-04 

EM ( 2 1) = -.741765E-04 DEL¢ 2, 1) = ~734907E-03 
EM( 2, 2) = -.853680E-04 

TRANSFER MATRIX PHI¢ 4) TRANSFER MATRIX DELTAC 4) 
EM( 1, 1) = ~235559E-04 DELC Ly 1) = ~118773E-04 
EM( 1, 2) = ~259254E-05 

EM( 2, 1) = -.259254E-05 DELC 2, 1) = -.164710E-04 
EM( on 2) 5 ~130299E-05 

TRANSFER MATRIX PHI 5) TRANSFER MATRIX DELTAC 5) 
EM( i, 1) = »748663E-06 DELC 1s. 1) = -344110E-06 
EM( Te 2) 5 ~6514565E-07 

EM( 2, 1) = -.651465E-07 DEL 2, 1) = .217079E-06 
EM (¢ 2, 2) = -.290989E~-07 
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TRANSFER MATRIX PHI ( 6) TRANSFER MATRIX DELTAC 6) 
EM( 1, 1) = ~197624E-07 DELC 1, 1) = -7359099E-08 
EM( 1, 2) 5 »150533E-08 
EM ( 2, 1) = -.150533E-08 DEL( 2, 1) = -.367553E-08 
EM ( 2, 2) = ~218458E-09 
GIVE THE INITIAL STATE (X(1,1)=---) 
x(1,1)=0.,0.* 
TIME x(1) X(2) 
+79 -OQOQ0E 00 -O000E 00 
1,57 -2929E 00 -7071E 00 
2.36 -1008E O01 -39692E 00 
3.14 -1758E 01 -5864E 00 
3.93 «2125E 01 -.2538E 00 
4.71 -1889E O01 -~.1100E 01 
5.50 «1158E 01 -.1478E 01 
6.28 -3207E 00 -.1159E 01 
7.07 -.1582E 00 -.2928E£ 00 
7.85 ~3256E~02 -6571E 00 
8.64 -7401E 00 -1163E 01 
9.42 »1663E 01 -9241F 90 
19.21 ~2263F O1 ~4467E-01 
11.00 -2192E O1 -.1009E 01 
11.78 -1462E 01 -.1654E 01 
12,57 ~4567E 00 ~.1514E 01 
13.35 ~.2735E 00 -.6351E 00 
14.14 -.3088E 00 +5194E 00 
14.92 »3980E 00 -1315€ 01 
15.71 -1481E 01 -1292E O1 
16.49 -2349E O1 -4317E 00 
17.28 »2503E O1 -.8138E 00 
18.06 ~1842E 01 -.1775E 01 
18.85 ~6889E 00 -.1890EF O01 
19.63 ~.3244E 00 ~.1067E O1 
20.42 -.6256E 00 -2718E 900 


END OF EXECUTION 
TO CONTINUE, GO TO THE TOP OF A NEW PAGE 
AND PRINT AN ASTERISK 


Figure 5,18 Console transaction for example 5.3 


when time delay = 7 min 


U: unit step 


Figure 5.19 Response curves of system for example 5.3 when time delay = 


a 


4 


eZ 
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CHAPTER 6 


COMMENTS AND SUGGESTIONS FOR FUTURE RESEARCH 


In obtaining Fas by the use of a digital computer the virtues of 
the series expansion technique are its simplicity and ease in programming. 
It is not necessary to find the eigenvalues of A. There is, however, some 
computational disadvantage to the series expansion method, This comes from 
the convergence requirements for the series. In general, it is reasonable 
to compute a by the power series when T is small, The running time for 
the matrix exponential simulation will be among the longest of various 
schemes. Use of the Jordan Canonical form, for example, requires 
considerably more programming, but will run in a fraction of time needed 
for the series solution, 

Some suggestions concerning the bound on the error in the 
evaluation of the matrix exponential when the matrix A is known with some 
error are given by Levis (10), 

The simulation technique for linear time-invariant dynamic 
systems has been tested, and it was found that the use of the augmented A 
matrix (X(t) = A X(t) + D U(t) can be expressed as v(t) - [¢ p|uce where 
V(t) -§] ) greatly improved the procedure. The reason is that the actual 
reduction of the elements of the augmented matrix times T to values less 
than one can be performed successfully. However, this method cannot be 
used for calculating the digital version of the control transition matrix. 

Another scheme that can be used to check the error bound in the 
state is to divide the time region of interest in two or three parts. 
Preferably these times should be powers of two times the sampling time. 


Next, compute the matrix exponential at the desired sampling time. 
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Recursively multiply it until the matrix exponential is found for the other 
selected times. The state at those times can be found and saved. Now, 

using the recursive process of state evaluation at the sampling time, 
compare the state with the selected ones, If the error is unacceptable, 

the state with less error can be used as a new initial condition, and the 
procedure may be continued. 


It was found in chapter 3 that the elements C form an array 


1,j 
of infinite order, The first row is of main importance because its elements 
are the terms of a, Therefore, the truncation technique already discussed 
can be used. 


In a similar fashion, the elements are actually the terms of 


Cyt 
the infinite series -. It is reasonable to expect smaller values of these 
norms as "i" grows, Therefore, intuitively the number of terms used to 
truncate the first row can be used to truncate $,(T), $2(T), etc. It would 
be interesting to make a study about how the truncation terms should be 


taken in each row in order to save computation time while maintaining 


accuracy, 
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TRANS 


Purpose: to compute the time response of linear time-invariant 


systems. 


Inputs: order of system (M = ); sampling time (T = ); final 
time (TF = ); number of input signals (R = ); the 


augmented A matrix and the initial state (X(1) = ). 


Outputs: the transition matrix; the current time; and the state 


of the system. 


Remarks: main program. Subroutine called by TRANS: EXPMAT, and 


DISTUR. 


MAGDA 


LUPE 


ALICIA 


JULIA 


ALMA 


FANNY 


MARTA 


OLGA 
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PROGRAM COMMON Ay EMs Ms RIds Ro X 

DIMENSION X20} 9YC2U) 9E( 20) sPE( 20) oX1 (20) 
DIMENSION EMP (400 9H) 9AC 400 9H) 9EM( 400 9H) 
INTEGER I sJeMeRoWI SH 

FORMAT VARIABLE FM 

VECTOR VALUES H=29190 

PRINT COMMENT SGIVE ORDER OF SYSTEM (M = )S& 
PRINT COMMENT $SAMPLING TIME (T = }» FINAL TIME (TF = 1% 
READ DATA 

PRINT COMMENT & & 

FM=M 

PRINT COMMENT $IS THERE ANY DISTURBING SIGNALS 
RFAD FORMAT S3sWISH 

VECTOR VALUFS S3 = $ C3*S$ 

WHENEVER WISHeF eSYESS 

PRINT COMMENT $3 $ 

PRINT COMMENT SGIVE NUMBER OF INPUT SIGNALS (R = )$ 
READ DATA 

M=M4R 

OTHERWISE 

R= 

END OF CONDITIONAL 

H(2)=M 

PRINT COMMENT $ $ 

PRINT COMMENT SGIVE THE A MATRIX (Allo1)=--9A(2s1)}=--)4 
THROUGH LUPEs FOR I=1lslsIeGeM 

READ PATA 

PRINT COMMENT & & 

PRINT COMMENT S$GIVE INITIAL STATE (X(1)=--)5 
READ DATA 

THROUGH ALICIA» FOR I=1l5151eGe(M-R) 

XICT}=XC1) 

TA=T 

WHENEVER ReNFev 

EXECUTE DISTURe (TA) 

J=M-R+] 

THROUGH JULIAs FOR I=JsleleGeM 

XI(J)=X(J) 

CONT INUF 

END OF CONDITIONAL 

THROUGH ALMA» FOR I=lsloleGe(M) 

E(1}=06 

TZ=T 

EXECUTE FXPMATe(T) 

THROUGH FANNY» FOR I=lolsIeGeM 

THROUGH FANNY» FOR J=lslsJeGeM 

PRINT FORMAT CUATROsIsJsEM( I 9J) 

VECTOR VALUES CUATRO = SIH 9S893HEM(91491H991493H) =9E14e6*d 
CONT INUF 

THROUGH MARTAs FOR I=lslsleGe(M-R) 

THROUGH MARTAs FOR J=lslsJeGeM 

EMP( I sJ)=EMCI J) 

WHENEVER (M-R) el 6 

PRINT COMMFNT $ $ 

PRINT FORMAT Sls (T=lolsleGe(M~R)sI) 

VECTOR VALUES Sl = & »S6s4HTIMEsSBs'FM*(2HX(s1191H)9S12)/*$ 
END OF CONDITIONAL 

TRANSFER TO TERESA 

TA=TAFTZ 
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WHENFVER ReNF ed 
EXECUTE DISTURe (TA) 
FND OF CONDITIONAL 
TERESA THROUGH ELENAs FOR I1=lolsIleGe(M-R) 
PE(1)=C6 
Y¥(1) 206 
THROUGH MARTA» FOR J=lolsJeGeM 
YOLT)=YCLI+EMP CT» J)*X(5) 
PFC T)SCEMP( Ts J)4+RIU)*E (CS) 4RIS*X( II 4P ECT) 
MARTA CONT TNUF 
FLENA CONT INUF 
ENORM=0, 
THROUGH ROSAs FOR I=1sls1eGe(M-R) 
ROSA ENORM=ENORM+e ABS e (PFI) ) 
WHENEVER ENORMeGFell0eePe-07) 
T=TA 
EXECUTE EXPMAT«{T) 
THROUGH ROSANAs FOR I=lolsleGelM-R) 
PE(T)=06 
Y(T)=06 
THROUGH ESTHER»s FOR J=lslsJeGemM 
YOCTISYCTIFEMC Ts J)*XI (CU) 
PE(TI=PF(TI4RIU*XI( J) 


FSTHFR CONTINUE 
ROSANA CONT INUF 
OTHERWISE 


TRANSFER TO SARA 
END OF CONDITIONAL 


SARA THROUGH CARMEN s FOR I=lslsleGe(M-R) 
X(T)=YCT) 
CARMEN F(T)=PECT) 
J=M=-R+1 
THROUGH LILIA» FOR I=J9lo1eGeM 
LILIA E(J)=06 


WHENEVER (M-R) el e6 

PRINT FORMAT S2sTAsX(1)eeeX(M-R) 

VECTOR VALUFS S2 = & 9S4sF6e2s'FMI(S39E14e6) *S 
OTHERWISE 

PRINT RFSULTS TA 

PRINT RESULTS X(1)eeeX(M-R} 

END OF CONDITIONAL 

WHENFVER TAeLe TFs TRANSFER TO OLGA 

PRINT COMMENT SEND OF EXECUTIONS 

PRINT COMMENT $TO CONTINUEs GO TO THE TOP OF A NEW PAGES 
PRINT COMMENT SAND PRINT AN ASTERISKS 

REA DATA 

TRANSFER TO MAGDA 

END OF PROGRAM 


YUTPOBE : 


Kemarks - 


EXPMAT 


compute the matrix exponential. 


saproutine called by TRANS, 
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ELENA 


DIANA 


OLGA 


SARA 


ALMA 
ESTHER 


YOLIS 


GLORIA 


ROSANA 


MARTA 


EXTERNAL FUNCTION (T) 

PROGRAM COMMON As EMs Ms RIJe Rs X 
DIMENSION A(4009H) 9EM(4005H) » TERM(4008H) sNTERM( 400 9H) 
DIMENSION X(20)9B(4U0sH) 

VECTOR VALUES H=2s1:s0 

INTEGER KsloJstoMolLLsY¥sQ 

ENTRY TO EXPMATe 

H(2)=M 

THROUGH ELENA» FOR I=lslsleGeM 
THROUGH ELENAs FOR J=lslsJeGem 
B(ITsJ)FA(T oJ) 

BCT sJ)SRC Is JI¥T 

AMIN=RA( 191) 

THROUGH DIANA» FOR I=25915J]eGe¥ 
WHENEVER RB(IsI)eleAMINs AMINEB( Is T) 
CONTINUE 

FAC=EXPe (AMIN) 

THROUGH OLGAs FOR I=1lslsIleGeM 
B(IsI)=B(Is1)-AMIN 

Y=eABSe8(191) 

THROUGH SARA» FOR I=1lsolsoIeGeM 
THROUGH SARA»s FOR J=lslsJeGem 
WHENFVER eABSe(Bll9J)) eGelY¥t00) s Y=eABSe( Bll oJ)) 
CONT INUF 

TAP=le 

YE=V¥+O056 

THROUGH ALMAs FOR Q=1919QeGeld 
TAP=2 o* TAP 

WHENEVER TAPeGEeYEs TRANSFER TO ESTHER 
CONT TINUE 

Y=TAP 

THROUGH YOLISs FOR I=lslsIeGeM 
THROUGH YOLISs FOR J=lslsJeGeM 

BUT s JERI Ts) /(Y +00) 
TERM(T9J)=BCI sJ) 

LL=0 

MAXH=O 6 

MAXV=06 

THROUGH MARIAs FOR I=lolsIleGeM 
SUMH=O06 

SUMV=O8¢ 

THROUGH ROSANAs FOR J=1lslsteGeM 
SUMH=SUMH+eABSeTERM( I 9J) 
SUMV=SUMV+eABSeTERM(UsI1) 

CONTINUE 

WHENEVER SUMH »Ge MAXH sMAXH=SUMH 
WHENEVER SUMVeGeMAXV »MAXV=SUMV 
CONTINUE 

NORM=MAXH 

WHENEVER MAXVeLl «NORM sNORM=MAXV 
WHENEVER LLeNEoUs TRANSFER TO DELIA 
SOLO=NORM 

K=2¢*NORM 

WHENEVER Kele2s K=2 

IN=K/2 

VECTOR VALUES CINCO = $1H »92HK=s14*S 
THROUGH SUSANAs FOR I=lslsIleGeM 
THROUGH SUSANAs FOR J=lslsJeGeM 
UNIT=06 


SUSANA + 
ISABEL 


EVA 


LILIA 


AURORA 


DELIA 


JULIA 


MAGUE 
MARTA 


OLIVIA 
ALICIA 


CARMEN 
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WHENEVER JeEols UNIT=le 

EM( IT sJ)=UNIT+B(I 9J) 

CONTINUE 

WHENEVER LL eGEeKs TRANSFER TO GLORIA 
LL=LL+1 

THROUGH LILIAs FOR L=lslslLeGeM 
THROUGH LILIAs FOR I=lsloleGeM 
NTERM(Lo1)=O6 

THROUGH EVAs FOR J=lslsJeGeM 
NTERM( Lo T)=NTERM(L31)4+B( Lod) *TERM( JI) 
CONTINUE 
EM(LeoI)=EM(LoT}4+NTERM(Lo TI) /(LL+1le) 
CONTINUE 

THROUGH AURORA» FOR I=lslsleGeM 
THROUGH AURORA» FOR J=1lslsJeGeM 
TERMUL sJI=ENTERM( Tod) /( LL +10) 
TRANSFER TO ISABEL 

EPS=SOLO/(K+2e) 
RIJ=NORM*SOLO/ ( (K+le)¥*(1e-EPS)) 
THROUGH JULIA» FOR I=lsloleGeM 
THROUGH JULIAs FOR J=lslsJeGeM 
WW=eABSe (EM( I 9J)*100eP 0-7) 
WHENEVER RIJeGeWW 

K=K+IN 

TRANSFER TO ISABEL 

OTHERWISE 

TRANSFER TO JULIA 

END OF CONDITIONAL 

CONT TINUE 

THROUGH ALICIAs FOR LL=loslsllb eGeQ 
THROUGH MARTA» FOR L=lslsLeGeM 
THROUGH MARTAs FOR I=1lolsIleGeM 
TERM(LoI)=QO. 

THROUGH MAGUF» FOR J=l91lsJeGeM 
TERM(LsI)=TERM(L oI) +EM(L oJ) EMC JIT) 
CONT INUF 

CONTINUF 

THROUGH OLIVIAs FOR I=lelsIeGeM 
THROUGH OLIVIAs FOR J¥lslsJeGeM 
EMCI sJ)=TERM( I 9J) 

CONTINUE 

PRINT COMMENT $ & 

PRINT COMMENT $ TERMS OF THE MATRIX EXPONENTIALS 
THROUGH CARMENs FOR I=lslsIeGeM 
THROUGH CARMENs FOR J=lslsJeGeM 
EMCI sJ)=FAC#EM( I 5J) 

FUNCTION RETURN 

END OF FUNCTION 


Me 
Mery 


* 


+ 
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TIMDEL 


Purpose: to compute the time response of linear systems with 


lumped parameters and time delays. 


Inputs: order of system (M = ); sampling time (T = ); time 
delay (TD = ); final time (TF = ); number of input 
signals (R= ); the A matrix; the B matrix; the Dy 


matrix; the D, matrix; the initial state (X(1,1) = ). 


2 
Outputs: the plant transition matrices, the control transition 
matrices if desired; the current time; and the state 


of the system, 


Remarks: main program. Subroutines called by TIMDEL: DELFOR, 


and PERTUR. 


MAGDA 


MELA 


MALENA 


ALMA 


BERTA 
JULIA 


PROGRAM COMMON EMsDELF »MsRsWeAeBeD1 sD2 oU 

DIMENSION FM(40005H) sDELF(4C005H) sX(4009G) sA(4009G) 
DIMENSION B(4U05G) 921(4003G) sD2(4009G) sUI4U0sE) 
INTEGER LIsJeKobLoMoNoLl oZowoRsREL »MMoWISHs JJ 

FORMAT VARIABLE FM 

VECTOR VALUFS G=29190 

VECTOR VALUES F=29199 

VECTOR VALUES H#=3519030 

PRINT COMMENT S$GIVE ORDER OF SYSTEM (M = )& 


PRINT COMMENT SDESIRED SAMPLING TIME (T = 18 

PRINT COMMENT STIME DELAY (TD = )» FINAL TIMF (TF = )% 
RFAD DATA 

FM=M 


PRINT COMMENT & & 

PRINT COMMENT $I1S THERE ANY DISTURBING SIGNALS 
READ FORMAT S39WISH 

VECTOR VALUFS S3 = $ C3*$ 

WHENEVER WISHeEe SYESS 

PRINT COMMENT & $& 

PRINT COMMENT $SNUMBER OF INPUT SIGNALS (R = )$ 
RFAD DATA 

OTHERWISE 

R=0 

END OF CONDITIONAL 

RELSTD/T4+002 

Gt2)=M 

H(2)=M 

H(3)=M 

PRINT COMMENT & & 

PRINT COMMENT SGIVE THE A MATRIX (A(191)#-=-sAl291)=-~)$ 
THROUGH MELAs FOR I[=lsloTeGeM 

READ DATA 

PRINT COMMENT © $ 

PRINT COMMENT $GIVE THE B MATRIX (B(191)=--»9Bl291)=--)$ 
THROUGH MALFNA® FOR I=1l91ls1leGeM 

RFAD DATA 

WHENEVER ReFeOsTRANSFER TO JULIA 

PRINT COMMENT £ & 


PRINT COMMENT $GIVE THE D1 MATRIX (D1(191)=--sD1(2sl)5--)$ 
THROUGH ALMAs FOR I=1lslsIeGeM 

READ DATA 

PRINT COMMENT $ & 

PRINT COMMENT SGIVE THE D2 MATRIX (D2(191)=--9D2(291)=--)% 
THROUGH BERTAs FOR I=lslsIeGeM 

RFAD DATA 


EXECUTE DELFORe (T) 

PRINT COMMENT $DO YOU WISH TO HAVF THE TRANSITION MATRICESS 
READ FORMAT S3sWISH 

WHENEVFR WISHeFeSYFSS 

THROUGH DULCE» FOR L=lsleleGewW 

LL=L~1 

WHENEVER ReEoO 

PRINT FORMAT OCHOsLL 

VECTOR VALUES OCHO = $1H +S8s15HTRANSFER MATRIXsS29 
14HPHI(91491H) *$ 

OTHERWISE 

PRINT FORMAT SEISsLLUsLt 

VECTOR VALUES SEIS= $1H »S8215HTRANSFER MATRIXsS29 
14HPHI (91491H) sS8s22HTRANSFER MATRIX DELTA(s14s1H)*S 
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CFLIA 
GLORIA 


ALICIA 


MARTA 


SILVIA 
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CONT INUF 

WHENEVFR ZeFels L=L-1 

THROUGH ALICIAs FOR K#l91loKeGeW*REL 
THROUGH ALICIAs FOR I=lelelTeGeM 
BCKel)=X(K+t91) 

WHENEVER KeEeWs TRANSFER TO ALICIA 
A(KsT)2U(K+19T1) 

CONT INUF 

THROUGH MARTA» FOR K=1l91skeGeW*REL 
THROUGH MARTA»s FOR I=lolsleGeM 
X(KeT)=B(Kel) 

WHENEVER KeEeWs TRANSFER TO MAF 7* 
UKs TVSA(KoT) 

CONTINUE 

TRANSFER TO SONIA 

PRINT COMMENT SEND OF EXECUTIONS 
PRINT COMMENT STO CONTINUEs GO TO THE TOP OF A NEW PAGES 
PRINT COMMENT SAND PRINT AN ASTERISKS 
READ DATA 

TRANSFER TO MAGDA 

END OF PROGRAM 
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DELFOR 


Purpose: to compute the plant transition matrices and the 


control transition matrices. 


Remarks: subroutine called by TIMDEL. 


DELIA 
YOLIS 


MAGUE 


ROSANA 


MARTA 


SALOME 
ISABEL. 
FANNY 
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EXTERNAL FUNCTION (T) 

PROGRAM COMMON EMsDELF »MoRoWoAsBsD1sD2 9U 
DIMFNSTON €C(110009H) sA(4009G) 9B(4005G) sEM( 4000 ott) *XX(400sG) 
DIMENSION TERM(400sG) sNTERM( 400 9G) »UU(4009G) #91(400%G) 
DIMENSION. DELF (4000 5H) »D2(4009G) sU(400sE) 
INTEGER IsJeKoLoMoNsY¥sQeRoW 

VECTOR VALUES H=3919030 

VECTOR VALUES G=2s1:0 

VECTOR VALUES E=291%0 

ENTRY TO DELFORe 

G(2)=M 

H(2) eM 

H(3)=M 

LINDA=06 

ROSA==1le 

THROUGH YOLIS» FOR I=lolsIeGeM 

THROUGH YOLIS»s FOR J=lolsJeGeM 

WHENEVER JeGeRe TRANSFER TO DELIA 

DIC IsJISDI (1s J) *T 

DAC eJ)=D2C Ts JI*T 

ACT sJIHACT 9 J)*T 

TERM( IT sJISACT oJ) 

BOTs J eB Tsu) *T 

N=0 

MAXH=06 

MAXV=06¢ 

THROUGH MARIA® FOR I=lolsleGeM 

SUMH=0 6 

SUMV=06 

THROUGH ROSANAs FOR J=lslsJeoGeM 
SUMH=SUMH+eABSe TF ERM( 19 J) 
SUMV=SUMV+eABSeTERM(Js1) 

CONTINUE 

WHENEVER SUMH eGeMAXH »MAXH=SUMH 

WHENEVER SUMVeGeMAXV »MAXV=SUMV 

CONTINUE 

NORM=MAXH 

WHENEVER MAXVel eNORM sNORM=MAXV 

WHENEVER LINDAeNEeO0es TRANSFER TO CARMEN 
WHENEVER NeNE«Os TRANSFER TO CHELA 
SOLO=NORM 

K=2e*NORM 

WHENEVER KeLe2» K=2 

W=1 

IN=K/2 

THROUGH SALOMEs FOR I=1s1lsIeGeM 

THROUGH SALOMEs FOR J=1s1lsJeGeM 
CllslsJ)=06 

WHENEVER JeoEolIsCllsIsJ)=le 
EM(WsTeJ)#C(1lsleJ) 

XX(TsJ) 2EM(WelsJ) 

UUCTst)=06 

TERM(TsJI=C(lslTsJ) 

N=0 

WHENEVER NeGEeKeANDeLINDAcEsO0es TRANSFER TO MAGUE 
WHENEVER NeGEeKeANDeLINDAeNEsOes TRANSFER TO ELENA 
N=N+1 

THROUGH HILDA» FOR L=ls91lslLeGeM 

THROUGH HILDA» FOR I=1lselsleGeM 
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LILA 


HILDA 


JULIA 


CHELA 


EVA 
FLENA 


AURORA 
MARTA 


IRMA 
PATY 


SONTA 


CARMEN 


JOSEFA 
YoCO 


ELISA 


NTERM(LoI)=06 
THROUGH LILIAs FOR J=lolsJeGeMm 
WHENEVER LINDAck eOesCI(NtlsJol )=06 


NTERM( Ls L)=NTERM(L9 IT) +A(L oJ) #TERM( Je 1) +B(L oJ) *CINt 1 oT) 


CONTINUE 
EM(WeLoT)=EM(WeL sl) +NTERM(L9 IT) /(NtLINDA) 
WHENEVER ReFeOs TRANSFER TO HILDA 
XX (Lo DIEXXCLs DI4NTERM( LoD) 70 CNtL INDA) #UN+LINDA+16)) 
CONTINUE 

THROUGH JULIAs FOR I=1lso1lsleGeM 
THROUGH JULIAs FOR J=l»lsJeGeM 
C(N+l ols J)=NTERM( IT 9J)/(N+LINDA? 
TERM( Ts J)=C(N+1 91 9J) 

TRANSFER TO FANNY 

EPS=SOLO/(K+2e) 
RIJ=NORM*SOLO/((K+le)*(1e-EPS)) 
THROUGH EVAs FOR I=lolsleGeM 
THROUGH EVA» FOR J=lslsJeGeM 
WW=eABSe (EM( Wel oJ) *1000eP 0-07) 
WHENFVER RI JeGeWW 

K=K+IN 

TRANSFER TO FANNY 

OTHERWISE 

TRANSFER TO EVA 

END OF CONDITIONAL 

CONT INUF 

LINDA=LINDA+41l 6 

ROSA=ROSA+1le 

WHENEVER ReEeOs TRANSFER TO PATY 
THROUGH MARTAs FOR L=lslslLeGeM 
THROUGH MARTAs FOR I=lsloleGeR 
TERM(LoI)=00 

THROUGH AURORA» FOR J=1s1sJeGeM 
TERM(LoT)=TERM(L oI) +XX( Lod) ¥DL (Jo +HUU(L J) #D20 eT) 
CONTINUE 

CONTINUE 

THROUGH 'IRMAs FOR I=1lsls1leGeM 
THROUGH IRMA» FOR J=lslsJeGeR 
DELF (Wel sJ)=TERM( IJ) 

THROUGH SONIAs FOR I=lsloTeGeMm 
THROUGH SONIAs FOR J=191l9JeGeM 
TERM(IsJ)=EM(Wol oJ? 

UUCT sJI=XX¢I9J) 

TRANSFER TO MAGUE 

WHENEVER NORMeLEel0eePe-O7 TRANSFER TO DIANA 
THROUGH YOCOs FOR L=1lslsl«GeM 
THROUGH YOCOs FOR I=lolsleGeM 
NTERM(LoI)=06 

THROUGH JOSFFAs FOR J=1l2lsJeGeM 
NTERM(LoI)=NTERM(LoT)4+BlbLsJ)¥*C{(lsdel) 
CONT INUF 

W=W+1 

THROUGH ELISA» FOR I=1lslsIeGeM 
THROUGH ELISAs FOR J=1sleJeGeM 
ClLeTsJV=NTERMUT 9J)/(ROSA+1 6) 
XXOTetJ=COlsIsJ)/CROSA+2 6) 
EM(WsTsJ)=C(lsIsJ) 
TERM( I sJ)=CC1lsIsJ) 

TRANSFER TO ISABEL 


DTANA 


vee 
1 

t : 
TON 
Se oF 
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PERTUR 


Purpose: to compute the forcing signal vector at the current 


time, The program keeps track of the past, 


Remarks: subroutine called by TIMDEL. 


EXTERNAL FUNCTION (TAsLL) 

PROGRAM COMMON’ EMsDELF sMeoRsWsAsBsD1l2D29U 

DIMENSION EM(4000sH) sDELF(40009H) 9A (400 9G) 9B (4009G) 
DIMENSION 01(4009G)9D2(4009G) sU(4Q005E) 

INTEGER IslLL»sRoWoM 

VECTOR VALUES G=29190 


‘VECTOR VALUES E=29190 


VECTOR. VALUES H=391s090 
ENTRY TO PERTURe 

G(2)=M 

H(2)=M 

E(2)=W 

H(3)=M 

U(LL 91) er 
U(LL 92) =-+<- 


U(LL 9R ) =- == 
FUNCTION RETURN 
END OF FUNCTION 
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